Steady entanglement out of thermal equilibrium 
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We study two two-level atomic quantum systems (qubits) placed close to a body held at a tem- 
perature different from that of the surrounding walls. While at thermal equilibrium the two-qubit 
dynamics is characterized by not entangled steady thermal states, we show that absence of thermal 
equilibrium may bring to the generation of entangled steady states. Remarkably, this entanglement 
emerges from the two-qubit dissipative dynamic itself, without any further external action on the 
two qubits, suggesting a new protocol to produce and protect entanglement which is intrinsically 
robust to environmental effects. 

PACS numbers: 03.65.Yz, 03.67.Bg, 03.67.Pp 



Introduction. — Entanglement represents one of the key 
features in quantum mechanics pQ, both from a funda- 
mental point of view, due to its connection to non local- 
ity [2], and from an applicative one, due to its crucial 
role in quantum information [3]. Environmental noise |4] 
induces decoherence effects [5 and is typically responsi- 
ble for the fragility of entanglement |6] . This represents 
one of the major obstacles to the concrete realization 
of quantum technologies related to quantum information 
processing [TJ [3] . A huge effort has been dedicated to the 
comprehension of the detrimental environmental effects 
[6HT0] and in conceiving suitable approaches to contrast 
the natural decay of quantum correlations [11 . Several 
approaches have been theoretically proposed and experi- 
mentally tested. They include reservoir engineering [IT], 
feedback methods [12], distillation protocols [13], deco- 
herence free-subspaces [14], non-Markovian effects [7], 
weak measurements [15], quantum Zeno effect [16 , dy- 
namical decoupling [17] and reservoir monitoring [18] . 

Here, we introduce a new direct procedure to protect 
entanglement realized by bringing the environment of a 
two-qubit system out of thermal equilibrium. Physical 
systems consisting of two qubits in a common environ- 
ment in absence [19] or presence [20] of matter have been 
largely investigated at thermal equilibrium, pointing out 
the creation of entanglement due to the field mediated in- 
teraction, which however typically washes off asymptoti- 
cally. Efforts have been also done considering two qubits 
interacting with independent ideal blackbody reservoirs 
at different temperatures [21 . New possibilities emerg- 
ing in realistic systems out of thermal equilibrium have 
been recently pointed out in different contexts ranging 
form heat transfer [22 , to Casimir-Lifshits forces [23H25] 
and atomic dynamics [26] . 

In this Letter, we consider a system made by two qubits 
interacting with the complex electromagnetic field out of 
thermal equilibrium resulting from the presence of bodies 
at different temperature and whose scattering properties 
play a key role. We will show that this particular environ- 
mental noise has two remarkable effects: it contrasts the 
usual dechoerence between the qubits, and it generates 



steady entangled states. This production and protection 
of entanglement is obtained without any further exter- 
nal actions on the two qubits, such as the use of lasers 
or complex procedures involving measurements on the 
qubits or on the environment, and without the need of 
initializing the total system in a given configuration. 

Physical system and model. — We consider two qubits 
q= 1,2, whose ground \g) and excited \e) internal lev- 
els are separated by the frequency uo = u\—w x g =uj^—Ug, 
interacting with a complex environment consisting in a 
stationary out of thermal equilibrium electromagnetic 
field. This is the result of the field emitted by a body 
of arbitrary geometry and dielectric permittivity, held 
at the temperature Tm, and of the field emitted by far 
surrounding walls held at temperature Tw, eventually re- 
flected and transmitted by the body (see Fig. [llwhen the 
body is a slab) . The walls have an irregular shape and are 
distant enough from the qubits such that the field they 
would produce at the qubits position (in the absence of 
the body) would be a blackbody radiation independent 
from the walls composition and geometry [23] . The total 
Hamiltonian has the form H = Hs + He + Hi, where 



H s = Eg E^ 



>^>lm 



being a^, 



(n|, is the 



free two-qubit Hamiltonian and He the free environmen- 
tal Hamiltonian. The interaction between the qubits and 
the environment in the multipolar coupling and in dipole 
approximation is described by Hj = — Eo D g • E(R g ) 




FIG. 1: (color online). Two qubits close to a slab at temperature 
Tm different from the temperature of the surrounding walls, Tw- 



[27] , where D g is the electric-dipole operator of qubit q 
(being q (g\T) q \e) = d* 7 ), and E(R g ) is the electric field 
at the position H q of qubit q. 

Master equation. — The starting point to study the 
two-qubit dynamics is the von Neumann equation for the 
total density matrix, which in the interaction picture is 
Ptot(t) = -^[ifj(£),/0totCO]- B y tracing over the environ- 
mental degrees of freedom, after the Born, Markov and 
rotating wave approximations, the master equation for 
the reduced two-qubit density matrix becomes 



d P = AiHs + 8 s ,p]-iY / A^'MKKCp] 



dt 



h l 



+ J2 Tqq '("){°ip° q J - lwiK"K'eH,p}) (i) 

q,q' 
q,q' 

where 5s is an operator related to the level frequency 
shifts, not playing any role in the following. Function 
A qq (lj) represents temperature independent induced co- 
herent (dipole-dipole) interaction between the qubits, 
while T^ (±o;) are individual (q = q') and common field 
mediated collective (q =fi q') qubit transition rates, re- 
lated to both quantum and thermal fluctuations of the 
electromagnetic field at the qubits positions. 

To discuss the properties of master equation (IT]), we 
will use two standard different basis: the decoupled bases 
{|1> = \gg),\2) = \eg),\3) = \ge),\4) = \ee)}, and the 
coupled bases {|G> = |1),|A) = (|2) - |3))/V2, |S) = 
(|2) + |3))/>/2, |E) = |4)}, where the collective anti- 
symmetrical and symmetrical states |A) and |S) are com- 
binations of the decoupled states |2) and |3). 

X states and concurrence. — In the decoupled basis, 
master equation (fil implies that the dynamics of the ele- 
ments along the two main diagonals of the two-qubit den- 
sity matrix (forming an X-structure) is independent from 
that of the remaining ones. Then, an initial state with an 
X-structure maintains its form in time. Moreover, terms 
outside the two main diagonals, are washed off asymp- 
totically. Bell, Werner and Bell diagonal states belong to 
the class of X states [30], which arise in a wide variety 
of physical situations and are experimentally achievable 
[31]. In the following we will deal with X states. 

We quantify the two-qubit entanglement by means of 
the concurrence C(t) (C = for separable states, C = 1 
for maximally entangled states) [32]. For X states, using 
Pij = (*Mj)j ^ takes the simple form [33 



C(t) = 2max{0,K 1 (t),K 2 (t)}, 



Kt(t) = \p 23 (t)\ - y/pii(t) P4A (t), 



(2) 



K 2 (t) = |p U (t)| - Vp22(t) P33 (t). 



Eq. (fil) induces an exponential decay for pi4(t), so that 
in the steady state only K\(t) could lead to C(oo) > 0. 



Thermal equilibrium. — For Tw = Tm master equation 
([!]) describes the qubits thermalization towards the diag- 
onal thermal equilibrium state [34 



/pn(oo) \ 

P22(oo) 
P33(oo) 



^eq 



/ [l + n( W ,T)] 2 \ 
n(w,T)[l + n(w,T)] 
n(w,T)[l + n(w,T)} 

V n(co,Tf J 



(3) 



(e fc B^ 



where Z eq = [1 + 2 n(u;,T)] 2 and n(u,T) 
l) -1 . This state is universal, it depends only on the 
ratio hw/k^T, remaining insensible to all system details. 
Being \p 23 (oo)\ = 0, Ki(oo) is always negative, resulting 
in not entangled steady states. In terms of the density 
matrix in the coupled bases, using px = Pxx = (X|/o|X), 



\P23\ 



(ps - Pa) 2 + (pas - Pis) 2 



(4) 



which is equal to zero in the steady state since Pas(oo) = 
and ps(oo) = pa(oo). The latter identity has not to be 
valid out of thermal equilibrium, allowing Ki(oo) > 
hence producing steady entanglement. 

Out of thermal equilibrium. — For Tw ^ Tm, the anal- 
ysis of Eqs. ([l][2| is much more rich and delicate. The 
T and A functions depend on the correlation functions 
of the electromagnetic field, which at thermal equilib- 
rium can be directly evaluated exploiting the fluctuation- 
dissipation theorem (FDT). Out equilibrium, the FTD is 
not valid in general. Nevertheless, we assume that the ra- 
diation emission by the body and the walls has the same 
characteristics it would have at thermal equilibrium at 
the source temperature [23, 25]. This allows to compute 
the correlation functions by indirectly using the FDT, as 
recently used to study the dynamics of a single atom [26] . 
The transition rates in Eq. (fiT) can be set under the form 



T<* M ^V r o Wm{[1 + n(u>,T w )]a#(w) 
[l + n(w,T M )]a$'(w)} 



(5) 



r™'(-«) =VrgMrg'M{n( w ,Tw)a$'(u;r 

+ n(w,r M )]a«'M*}, 

where a«'(w) = Em^WV^MW, < M = 

£ M 4d\*[d 9 ']*«V)W, being [d% = [d«]i/\d*\, and 
To(cj) = |d g | 2 o; 3 /3fi7reoc 3 is the vacuum spontaneous- 
emission rate of qubit q. [a,y[ (<jj)]u> and [a^ (oo)]w are 
temperature independent functions, which depend on all 
the other system parameters (qubits positions and fre- 
quency, dielectric and geometric properties of the body) 
and can be expressed in terms of the reflection and trans- 
mission operators associated to the body (see appendix). 




FIG. 2: (color online). Scheme of the rate equations of Eq. (mi: 

T l(A) = r s(A)(i + "s(A)), r™ A) =r s(A) n s(A) . 



The A" (w) function of Eq. <Q can be expressed as 



A qq (oj) 



^/rg( W )rg'(a;) 



xP 



/ 



+00 , ,/3 



a Aia/ a% (a/) +ag?(y) 

27T CJ — C*/ 



(6) 



and calculated using the Kramers-Kronig relations (see 
Appendix). 

Analytical investigation. — To illustrate the new qual- 
itative and quantitative behaviour of the entanglement 
out of thermal equilibrium, we first consider an instruc- 
tive case allowing a direct interpretation. Let us con- 
sider r n (±cj) = T 22 (±oj) = T(±oj) andT 12 ( 21 )(±u;) e M. 
These conditions are verified for identical qubits (d 1 = 
d 2 = d) in equivalent positions with respect to the body 
(in the case the body is a slab, z\ = z^) and with d 
real and directed or along the z axis or along the x — y 
plane. In this case, master equation ([I]) implies in the 
coupled basis a set of rate equations for the populations, 
decoupled from the other density matrix elements: 



PG=r A (i + n A )PA + r s (i 

- (r A nA + r s ns)pG+, 
p A =r A n A p G + T A (1 + n A )p E 

ps =r s n s pG + r s (i + n s )pE - r s (i 

Pe =r A n A p A + T s n s ps 

-[r A (l + 2n A )+r s (l + 2n s )]p E . 



^s)ps 

r A (l + 2n A )p A , 
2n s )ps, 



(7) 



Here the derivates are with respect to To(oo)t [^0(00) 



.(i)/ 



^(2) 



Iq (oj) = Tq (oj)] , and we introduced the symmetric and 
anti-symmetric rates and the effective number of photons 



T A =a w (oj) - o<w(oj) + a M (u) - Om(u) 
Ts =a w (oj) + a^(oj) + a M (u) + «m M 



n A 



r A i 



%M - «wM] n (^> ^w) 



n s 



+ [a M (u) -«mM] n(cj,T M )| 

= jr{[ a w(w) +QwH] n(cj,T w ) 



(8) 



being a W (M)M = «w(M)(^) = a w(M)M- We remark 
that function A does not enter in the rate equations (|7|), 
which are depicted in Fig. [2] To each decay channel from 
|E) to |G), passing respectively trough |S) and |A), one 
can associate an effective number of photons ng( A ) con- 
fined between n(oj, Tw) and n(oj, Tm), which is equivalent 
to associate an effective temperature T S ( A ) confined be- 
tween Tw and Tm [26]. It is possible to show that while 
the coherences along the second diagonal decay exponen- 
tially to zero, the stationary solution of Eq. ffil is 



/ Pg(oo) \ 

Pa(oo) 
ps(oo) 
V Pe(oo) 



v neq 



neq 



/ (l + n A ) 2 (l + 2n s )r A + (l 

n A (l + n A )(l + 2n s )r A + [n A (l - 
n s (l + n s )(l + 2n A )r s + [n s (l + 
n A (l + 2n s )r A + (l 



V 



f2n A )(l + n s ) 2 r s 
-2n s )+n|(l + 2n A )]r s 
2n A ) + n A (l + 2n s )]r A 
f 2n A )n|r s 



(9) 



where Z neq is the sum of the elements of the vector on 
the right side of the above equation. Equation (|9J) , which 
reduces to (p3| for Tw = Tm, shows that out of equi- 
librium it is possible that ps(oo) ^ p A (oo), and implies 
|P2s(oo)| = \n s - n A \(T s + r A )/2Z neq trough Eq. @. 
This leads to the possibility to have Ki(oo) > in Eq. 
( 2 ) , corresponding to stationary entanglement. Using Eq. 
(9) in Eq. pi), we obtain for the steady concurrence 
C(oo) = 2 max{0,ifi(oc)}, with 



tfi(oo) = 



|n s -n A |(r s + r A )/2 



V(l + n A ) 2 (l + 2n s )r A + (1 + 2n A )(l + n s ) 2 r s 



ni(l + 2n s )r A + (l + 2n A )n 2 r s 



(10) 



which tends to zero at thermal equilibrium when ns = 
n A . Simplifying Ts, C(oo) becomes function of only 




n s 0.15 



FIG. 3: (color online). Part (a): steady concurrence [C = C(oo)] 
vs ns and tia for a fixed value Ta^s ~ 2.8 x 10 -4 . Part (b): 
maximum of concurrence, C m ax as function of Ta/Ps- 



Ta/Fs, ns and n\. We discuss this dependence in Fig. 
[3] (a), where C(oo) is depicted as a function of ns and 
nA for Ta/Fs « 2.8 x 10 -4 . It is shown that large val- 
ues of steady concurrence are obtained when the num- 
ber of photons associated to the two decay channels (see 
Figl2|, that is their effective temperatures, are enough 
distant between them. This physically corresponds to 
largely populate the antisymmetric state with respect to 
the symmetric one [see Eqs. (J9J) and Q]. By increasing 
too much nA at fixed ns, the steady entanglement starts 
to decrease (not shown in the figure). Part (b) shows 
that the maximum value of C(oo) reachable by varying 
ns and nA at a fixed value of Ta/Fs is 1/3, obtainable 
in the two cases Ta/Fs — » or Ts/Fa — > 0. We remark 
that up to now our findings do not rely on the specific 
choice of body's geometry or dielectric properties. 

Numerical investigation. — In order to discuss the prop- 
erties of C(oo) besides the case studied above, we solve 
Eq. ([!]) for the case where the body close to the two 
qubits is a slab of thickness <5, as depicted in Fig. [I] In 
this case a^ and a™ of Eq. ([5| can be expressed as in- 
tegrals over propagative and evanescent sectors (see Ap- 
pendix). We also choose a SiC slab, describing its dielec- 
tric permittivity with a Drude-Lorentz model, with a res- 
onance at uj r = 1.495 x 10 14 rad s _1 and a surface phonon- 
polariton resonance at u p = 1.787 x 10 14 rads _1 . Hence, 
relevant length and temperature scales are c/uo r ~ 2um 
and Huj r /kB — 1140 K. In Fig. [41(a) we plot concurrence 
of Eq. pi as a function of z 2 and Tm in the case of two 
identical qubits having electric dipole perpendicular to 
the slab, for fixed values of z\ = l/xm and Tw = 30 K. 
The plot evidences a large zone in the space of the pa- 
rameters corresponding to the generation of steady en- 
tangled states. The maximum value of C(oo), obtained 
for z\ i- z 2 = 1.28 and T M « 1300 K, is « 0.224. The 
characteristic time to reach this entangled steady state 



is 



io 3 [r (a 



[see part (b)]. The white line corre- 



sponds to the case z 2 = z\ and hence can be described 
by Eq. ( p~QJ ) . The maximum along this curve, obtained for 
Tm « 1200 K, corresponds to the red point in Fig. [31(a), 
being n A « 1.53 (T A « 680 K) and n s « 0.02 (T s « 90 
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FIG. 4: (color online). 
vs Z2 and Tm ■ Here z\ 



Part (a): steady concurrence [C = C(oo)] 
ljum, r\2 = 0.25/im (the qubits are distant 



V\ 2 + {zi - z 2 ) 2 } 1/2 ), T w = 30 K, S = 0.01/im, and u = 0.3ou r 
Part (b): dynamics of concurrence as a function of Fo(uj)t for three 
different initial conditions, the antisymmetric (green dashed line), 
the symmetric (blue dotdashed line) and the thermal state at 30 K 
(red solid line). 



K). The relevant difference between Ts and Ta is respon- 
sible of the high value of concurrence, « 0.217 (see also 
Fig. [3j) . However, by further increasing Tm (comport- 
ing an increase of the difference between Ts and Ta) the 
concurrence decreases. In Fig. El (b), by using the pa- 
rameters corresponding to the maximum of Fig. RJa), we 
show the time evolution of concurrence for three differ- 
ent initial states: the maximally entangled antisymmet- 
ric and symmetric states, and the not entangled thermal 
state at T = 30 K. We note how the protection of en- 
tanglement and its steady production, respectively, are 
independent on the initial two-qubit state. A systematic 
study shows also that by increasing the slab thickness (5, 
or the value of 7*12, or moving the atomic frequency uj 
towards the slab resonances, or changing the two qubits 
electric dipole orientations steady entanglement typically 
reduces. We observe that even a small amount of mixed 
state entanglement, here produced, could be then dis- 
tilled into a pure entangled state [35] . 

Conclusions. — We investigated the dynamics of two 
qubits interacting with a common stationary field out 
of thermal equilibrium. We predicted the occurrence of 
steady entangled states not depending on the initial two- 
qubit state, consisting then in a creation and/or pro- 
tection of entanglement according to the nature of the 
initial configuration. For a relevant class of parameters 
we derived an analytical expression for concurrence, and 
explained the entanglement production in terms of rate 
equations driven by two different effective temperatures 
associated to the two decay channels governing the pas- 
sage from the two-qubit excited state to the ground state. 
We numerically studied the case where the body close to 



the qubits is a slab, finding concurrence up to « 1/4 for 
conservative parameter s's values. While at thermal equi- 
librium the entanglement decays to zero faster if the tem- 
perature is increased, this new strategy to create and/or 
protect entanglement can be realized, quite counterintu- 
itively, starting from a thermal equilibrium configuration 
and increasing only one of the two temperatures of the 
system. To further increase the amount of steady entan- 
glement, systematic studies exploiting different body's 
geometries are envisaged. 



Appendix: a functions 



The expressions for the functions [a^ (uj)] 



a M MW appearing in Eq. 



and 

51, in the case the body 
forming part of the environment of the two qubits has an 
arbitrary geometry and dielectric function, are [the two 
qubits are placed in (ri, z\) and (r2, £2), where 1*1 and v<i 
are vectors in the xy plane] 



M] 



3irc 



E 



d 2 k f d 2 k 



2w ^ J (2tt) 2 J (2tt) 2 



D *(k-r„ 



«- k '-V){p,k|{e i ( fc **«- fc * , v)[e+(k,a;)] i [eJ(k',a;)]t, 

(TP^T* + llV^rt) + e < (*"« +fc **v)[g+(k ) w)] i [v(k' ) w)]?,^^ 
■e-« k ''* +k *''<'^{k,w)]i[e${V,^ 



(11) 



K MW = 



3ttc 



Y [—[ 



d 2 k' 

(2lt) 2 



z(k-r q -k' 



'■V)( p ,k|{e < ( fc **«- fc « , v)[g+(k, w )] i [ej(k' ! 



u 



^(pw) _ nV (V^) n \ + ^H _ p(ew)^ t _ r:p (pw) r t^ j | p , } k ^ 



where the operators 7£ and T are the standard reflection 
and transmission scattering operators, explicitly defined 
for example in [25], associated in this case to the right 
side of the body. They connect any outgoing (reflected or 
transmitted) mode of the field to the entire set of incom- 
ing modes. In the two previous equations, each mode of 
the field is identified by the frequency cj, the transverse 
wave vector k = (k x , k y ), the polarization index p (taking 
the values p = 1, 2 corresponding to TE and TM polar- 
izations respectively), and the direction or propagation 
<j) — ±1 (shorthand notation (j) — ±) along the z axis. 
In this approach, the total wavevector takes the form 
K^ = (k, (ftkz), where the z component of the wavevec- 



tor k z is a dependent variable given by k z = \ r^r — k 2 , 
where k = |k|. For the polarization vectors appearing in 



Eq. (11) we adopt the following standard definitions 



e^ E (k,cj) = z x k= -(-k y k + k x y), 



c£ M (k,u;) = -€^ E (k,cj) x K 4 



-(—kz + <fik z k), 



(12) 



where x, y and z are the unit vectors along the three 



axes and k = k/fc. In Eq. (11) we have also used 



<p,k|p(p w / ew V,k') = ^<p,k|n(p w / ew y,k'), (13) 



I[( pw ) and II^ ew ^ being the projectors on the propaga- 
tive (ck < oj, corresponding to a real k z ) and evanescent 
(ck > uo, corresponding to a purely imaginary k z ) sectors 
respectively. 

In the case when the body is a slab, we have at dispo- 
sition simple expressions for 1Z and T. As a result of the 
translational invariance of a planar slab with respect to 
the xy plane, its reflection and transmission operators, 
1Z and T, are diagonal and given by 



(p,k\K\p\k') = (2n) 2 S(k-k')S pp/ p p (k,u J ), 

{p,k\T\p',k') = (27T) 2 S(k-k')S pp ,T p (k,0j), 



(14) 



where the Fresnel reflection and transmission coefficients 
modified by the finite thickness S are equal to 



1 - e 

p p (k,cj) = r p (k,cj) _^ 2 



ZlK zrn d 



T p (k,uj) 



p {k,Lo)e 2ik *™s- 

t p {k,uj)i p (k,uj)e i( ~ k ^- k ^ d 
l-r 2 (k,co)e 2ik ^ s ' 



(15) 



In these definitions we have introduced the z component 
of the K vector inside the medium, 



= \ e(u) 



OJ* 



■k2, 



(16) 



s(lj) being the dielectric permittivity of the slab, the or- 
dinary vacuum- medium Fresnel reflection coefficients 



^TE 



^TM 



s(uj)k z 



s(uj)k z + k z 



(17) 



as well as both the vacuum-medium (noted with t) and 
medium-vacuum (noted with t) transmission coefficients 



txE 



txE 



2k z 



2k z 



tTM = 
?TM = 



2y/^j)k z 

s(uj)k z + k zm '' 

2^e(uj)k zm 
e(u)k z + k zm ' 



(18) 



Using Eq. (14), Eq. (11) reduces in the simpler case 



when the body is a slab to 

[A«'( W )]J i , + [B«'(a;)]«/+2[C«'(a;)] i 



x w 



M]i 



K2 MW = 



[A™ (w)]«, - [B« (w)]«/ + 2[D« (w)]i 



(19) 



where we have introduced the integral matrices 



[A«{ U )] W ~J2JJ ^e tt -(««-v)[JvW (A . >w)] +- 



+ + 



[£ gg (")]* 



3c 



E 



kdk 



4uj ^— ' J k z 
p u 



e *fc.(^-V)[JV«'(ife,a;)]+ + 



x (|p p (fc,w)| 2 + |r p (A;,w)| 2 ), 

3c 

4w 



[C«'(o;)]«' =|^E /' ^Re{e ifc *(*«+V) 



o fc* 



x [JV«'(fc,a;)]+-p p (fc,a;)}, 



[£>« (w)]i 



3c 

~4w 



E 



+oo 



^L -Im(fc,)( z ,+ V ) 

Im(^) 



x[7V«(fc, W )]++Im[ /0p (fc, W )], 



(20) 



where [we choose the x axis along the vector r q — ry 
whose coordinates in the plane xy are then (r^/,0)] 



[Nf {k,uj)\if 



^Ji(kr qq ,) 0> 

j^Ji(kr qq ,)-2J 2 (kr qq ,) | , 

qq 0; 



[JV?* (*,")]# = 



( 2 tt C ,!!k l J i( fer V) ~ kr qql J 2 (kr qq >) 



\ -i(j) , '^-^J 1 {kr qq >) 



-i^-^J^kr^ 



kr qq ,u>* J nKr qqf ) 




2c z k A 



Jv{kr qq >) j 



where i runs over the rows and %' over the columns and 
where J n (x) is Bessel function of the n kind. 

Concerning A qq (uj) of Eq. (|6|, we can exploit the con- 
nection between a functions and the imaginary part of 
the Green function of the system, 



Im Gut (Rq , Rg' , id) 



y qq 



(co)] ie + [a™ M], 



37reoc 3 



(21) 

where z and z' refer to the cartesian components of the 
field and G^/(R g , R g /,o;) is the zz' component of the 
Green function of the system, solution of the differential 
equation (for two arbitrary R and R') 



Vr x V R 



uj 



e(o;,R) 



(R,R',cj) 



6qC< 



■I^R-R') 

(22) 



being I the identity dyad and e(cj, R) the dielectric func- 



tion. Using Eq. (21), Eq. (6) becomes 

AM 'H = -^E™<i 9 V 



xP r^ imG»,(R 9 ,iyx) (23) 

-^X)[ d *]<[d*']i'ReG«/(R„IV,w), 



where Kramers-Kronig relations connecting real and 
imaginary parts of the green function have been used 
to compute the principal value of the integral. 

In the case of a slab, it can be shown that previous 
equation reduces to 



A«'( W ) =Ag«'( w ) + VrgMrg'( w ) £[a*]?[d* ],, 



i,i' 



(24) 



([cf mw - [z>r (, 



where [d ]$ = [d 9 ]a/|d 9 | and Aq 7 (lj) is the free term re- 
maining in absence of matter around the two qubits 

a«'m = - 



rg( W )rg'(o;){[d«*.d«'-(d«*.r) 



(d'-r)] 



■2(d«* -f) (d« •?) 



(f 2 — 1) cos f — f sin f 






cos r + r sin r 



}■ 



(25) 



being f = |r| = |r g — iy|c<;/c, and where we have intro- 
duced the new integral matrices 



[cr Hh 



3c 



£ 



8cj ^^ ./o kz 



kdk . 



x [N^'(k,Lo)}±r Pp (k,co)), 



(26) 



-Im(fe z )(z q +^/) 



x[A^(fc, W )]++Re(p p (fc,a,)). 
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